SUBROUTINE grid(x,xmin,xmax,s,n)
! Purpose: Generate grid x on [xmin,xmax] using spacing parameter s set as follows:
! s=1		linear spacing
! s>1		left skewed grid spacing with power s
! 0<s<1		right skewed grid spacing with power s
! s<0		geometric spacing with distances changing by a factor -s^(1/(n-1)), (>1 grow, <1 shrink)
! s=-1		logarithmic spacing with distances changing by a factor (xmax-xmin+1)^(1/(n-1))
! s=0		logarithmic spacing with distances changing by a factor (xmax/xmin)^(1/(n-1)), only if xmax,xmin>0
	IMPLICIT NONE
    INTEGER, INTENT(IN) :: n,s
	DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: x(n)
	DOUBLE PRECISION, INTENT(IN) :: xmin,xmax
	DOUBLE PRECISION :: c ! growth rate of grid subintervals for logarithmic spacing
	INTEGER :: i
	FORALL(i=1:n) x(i)=(i-1)/real(n-1)
	IF (s>0) THEN
		x=x**s*(xmax-xmin)+xmin
		IF (s==1) THEN
!			PRINT '(a,i8,a,f6.3,a,f6.3,a)', 'Using ',n,' equally spaced grid points over domain [',xmin,',',xmax,']'
		ELSE
!			PRINT '(a,i8,a,f6.3,a,f6.3,a,f6.3,a)', 'Using ',n,' skewed spaced grid points with power ',s,' over domain [',xmin,',',xmax,']'
		END IF
	ELSE
		IF (s==-1) THEN
			c=xmax-xmin+1
!		ELSEIF (s==0.0D0) THEN
!			IF (xmin>0.0D0) THEN
!				c=xmax/xmin
!			ELSE
!				STOP 'grid: can not use logarithmic spacing for nonpositive values'
!			END IF
		ELSE
			c=-s
		END IF
!		PRINT '(a,i8,a,f6.3,a,f6.3,a,f6.3,a)', 'Using ',n,' logarithmically spaced grid points with growth rate ',c,' over domain [',xmin,',',xmax,']'
		x=((xmax-xmin)/(c-1))*(c**x)-((xmax-c*xmin)/(c-1))
	END IF
END SUBROUTINE grid

